Explosive synchronization enhanced by time- delayed coupling 
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This paper deals with the emergence of synchronization in scale-free networks by considering the 
Kuramoto model of coupled phase oscillators. The natural frequencies of oscillators are assumed 
to be correlated with their degrees and a time delay is included in the system. This assumption 
allows enhancing the explosive transition to reach a synchronous state. We provide an analytical 
treatment developed in a star graph which reproduces results obtained in scale- free networks. Our 
findings have important implications in understanding the synchronization of complex networks, 
since the time delay is present in most systems due to the finite speed of the signal transmission 
over a distance. 

PACS numbers: 05.45.Xt, 89.75.Fb, 02.10.Ox 



INTRODUCTION 

Synchronization is an emerging collective phenomenon 
resent in science, nature, social life and engineering 
When a large population of limit-cycle oscillators is 
coupled, with slight different natural frequencies, a syn- 
chronous state emerges and those oscillators run at an 
identical frequency [l|,y|. For instance, experimental stu- 
dies have verified that neurons in the central nervous sys- 
tem synchronize their oscillatory responses [HI . Similarly, 
the emergence of synchronization has also been observed 
among electrochemical oscillators [6]. Indeed, the onset 
of synchronization has been verified in circadian rhythm, 
communication networks, power grids, social interacti- 
ons, cortical networks and ecology 

Winfree proposed a mathematically tractable model to 
explain the emergence of collective synchronization 
Although his approach allows describing the synchroni- 
zation dynamics in a population of coupled oscillators, it 
assumes that every oscillator feels the same mean-field, 
which is hard to be observed in real-world complex sys- 
tems. A more suitable model was proposed by Kura- 
moto, who considered oscillators coupled by the sine of 
their phase differences and phase oscillators at arbitrary 
frequencies jl, [9|. Therefore, each oscillator i obeys an 
equation of motion given by 



dt 



N 

uji -\- X Aij sm{Oj 



1,...,7V, (1) 



where A is the coupling strength, uJi is the natural fre- 
quency of oscillator i (generally distributed according to 
some function g{uj)), and Aij are the elements of the adja- 
cency matrix A, which represents the topology of a com- 
plex network. More specifically, elements Aij = 1 if two 



nodes i and j are connected, and Aij = 0, otherwise |ld |. 
In the Kuramoto model, each node tries to oscillate in- 
dependently at its own frequency, while its direct neigh- 
bors tend to synchronize it to all the other oscillators. 
It has been verified that coupling strengths higher than 
a determined threshold Ac produce the onset of synch- 
ronization [4]. The level of synchronization of the whole 
system is measured by the macroscopic complex order 
parameter 



r{t)e 
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where < r{t) < 1 measures the phase coherence of 
populations. When r{t) = 1, the system reaches the 
complete synchronous state and when r{t) =0 it reaches 
an asynchronous state. 

The Kuramoto model has also been analyzed conside- 
ring some modifications, such as different types of cou- 
pling functions and frequency distributions (J |TH , in or- 
der to enhance the synchronization level. For instance, 
Gardenes et. al. (lH investigated the Kuramoto model in 
complex networks by defining the natural frequencies of 
each node i as degree /c^, ie., uJi = ki. In contrast to the 
canonical choice for the frequency distribution, they ob- 
served a first-order phase transition to the synchronous 
state in scale- free networks, showing that this effect is due 
to the positive correlation between the nodes frequencies 
and the network topology. 

The Kuramoto model has also been modified to in- 
clude time-delay {e.g. |Tl-[l3|)- Time delay is observed 
in communication systems due to the finite speed of the 
signal transmission over a distance. For instance, the 
transmission of signals through unmyelinated axons pre- 
sents an 80-ms delay for propagation through a corti- 
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Figura 1: Coherence diagram r(A) for BA model with N 
1000 and (/c) = 6. 



cal network [17]. As verified by Dhamala et al. |18j . 
the time delay can enhance the network synchronization. 
They studied neural synchrony by taking into account the 
Hindmarsh-Rose neurons and observed that an extended 
region of stable synchronous activity could be achieved 
with low coupling strengths only if the delay were inclu- 
ded in the system. The contribution of the time delay 
to network synchronization was also verified by Perez et 
al The authors incorporated it in the Kuramoto 

model and verified that it was possible to obtain a per- 
fect synchronization in directed networks by relating the 
topology of the network and the phases and frequencies 
of the oscillators. 

In the current paper, we took into account the analyse 
by Gardefies et. al. 12\ and Choi et al [l3| altogether. 
We considered the natural frequencies of each node i as 
the node degree uji = ki and included a time delay in 
the Kuramoto model. Therefore, the set of equations of 
motion for N coupled oscillators is given by 



dt 



N 



\Y,A,,sm[Oj{t-T)-0,(t)]. (3) 



Note that the phase of node i interact with the phase of 
node j with a time delay r. 

We analyzed this modified Kuramoto model (Eq. [3]) 
in Barabasi- Albert (BA) networks [ll|. Fig. [1] shows the 
synchronization diagram r(A) taking into account a scale- 
free topology with N = 10^ and {k) = 6. As performed 
in [lH we increased the value of A adiabatically and com- 
puted the stationary value of the global coherence r for 
each value Aq, Aq + (5A, . . . , Aq + n5\. Fig. [1] shows 
that the presence of a time delay in the phase interac- 
tions allows reaching the synchronous state with smaller 
couplings than without time delay. In addition, it is in- 
teresting to observe that the system presents different 



phase transitions for each time-delay r adopted. For ins- 
tance, for r = 1.0, the transition is of first order, while 
for r = 0.3, r increases almost linearly with A and no 
phase transition is verified. 

In order to locally investigate the effect of the time de- 
lay on BA networks, a star graph was considered as an 
approximation of the scale-free topology. A star graph is 
a special type of tree, i.e., a central node connected with 
K leaves. Thus, for a star network with = i^+1 nodes, 
the peripheral node has degree ki = 1 {i = 1, ...,i^r) and 
the central node has kn = K connections. By using this 
approximation we studied the influence of hubs (nodes 
with higher frequencies) on neighbors nodes, which are 
expected to have lower frequencies, since the average clus- 
tering coefficient tends to zero in BA networks [20|, l21 |. 
As in the work by Gardehes et. al. [lH, we considered the 
natural frequencies of each node equal to its respective 
degree. Setting the rotated frame with the average phase 
V^(t) = '0(0) + l^t, where ft is the average frequency of 
the system, we defined the new variables as (f)i = Oi — Qt 
for the peripheral nodes and (j)H = 9h — for the hub. 
These definitions allowed rewriting Eq. [3] into the set of 
equations given by 



K 



OH 



(cj^ - 1]) + A ^ sin (t-r)- c^h (t) - ^r] , (4) 



3 = 1 



(j)i = {uj- n) + Asin [(t)H{t - r) - (t)i{t) - Qr] , (5) 

where we set ^/^(O) = without loss of generality. The 
global parameter defined in Eq. [2] can be rewritten into 



K + 1 



(6) 



Multiplying both sides of Eq. [6]by e '^i^H-\-^r) taking 
the imaginary part, we can rewrite Eq. [Has 

= {^H -^)- Xr{K + 1) sm{(j)H + l^r) - A sin fir. (7) 

Imposing the phase locked solution <j)H = we obtained 

{ujh — ^) sin Qr 



sin(0// + Qr) 



(8) 



Note that when the time delay vanishes, r = 0, the same 
result verified in [1^ for the hub's phase is recovered 
{i.e. sin (j)H = {uh — 0)/\r{K + 1)). Considering Eq. [5] 
and solving it for cos in the locked regime 0^ = we 
obtained 



cos ( 



{yt — uo) sm.{(j)H — ^t) 
A 



sm 



Q.t)\ [A2 - (cj - 



(9) 



A 



This equation is valid only for A ^ |cc; — which leads 
to the critical coupling Ac = |cj — 1]| for the onset of 
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Figura 2: Coherence diagram r(A) for {a)K = 60, {h)K = 70, 
{c)K = 80 and {d)K = 90 
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Figura 3: Average frequency Q as a function of time delay r 
for {a)K = 60, {h)K = 70, {c)K = 80 and {d)K = 90 



synchronization. Nevertheless it is important to stress 
that this value for the critical coupling depends on the 
time delay r, once the average frequency also depends 
on it (ie., Q = ^{r)). When r = 0, the result Q = 
2K/{K + 1) verified in [12] is recovered. Moreover, it 
is straightforward to show that at the critical coupling 
A = Ac the phases of the peripheral nodes and the hub 
obey the following condition 

6^{r)=^H-^^ = ^^rQ{r). (10) 

The above equation shows that the difference between 
the phases of peripheral nodes depends directly on the 
time delay r, which means that it is possible to tune the 
time delay in order to obtain different phase couplings 
between the peripheral nodes and the hub. 

In order to obtain insights on how the presence of the 
time delay r is related to the onset of synchronization, 
we constructed the synchronization diagram for the star 
networks, as shown in Fig. O Comparing the results pre- 
sented in Figures[2]and[5jinset), we can see that the time- 
delay enhance the network synchronization. The values 
of K were selected by taking into account the expected 
degree for a hub in a scale- free network of N = 10^, as 
we considered the star graph an approximation of a BA 
network. Remarkably, Fig. [2] shows that star graphs pre- 
sent the same behavior observed in the scale-free network 
in the presence of a time delay (see Fig(T]). More speci- 
fically, depending on the value of r, the star networks 
undergo different synchronization transitions. It is inte- 
resting to note that, like in Yeung et. al. [l3|, a subcri- 
tical Hopf bifurcation of the incoherent state at r = 1 is 
observed. In addition, a supercritical Hopf bifurcation at 
r = 2 is verified in Figure [TJ for the BA networks. 



This fact is due to the dependence of the mean fre- 
quency ft on r, i.e., Q = rt{r). We analyzed such de- 
pendence in order to understand how the critical cou- 
pling Ac varies for different time delays. Fig. [3] shows 
the mean frequency ^7 as a function of r. It is in- 
teresting to note that the particular choice of the fre- 
quency distribution g{uj) = P{k) leads to a dependence 
of ^{r) that is different than that verified in (l^, in 
which the authors found that frequency ft depends in- 
versely on the time delay r for the Kuramoto model with 
all-to-all connectivity. Fig. [3] shows an oscillatory de- 
pendence between Q and r, which explains the different 
behavior of the synchronization diagrams displayed by 
the scale-free and stars networks in contrast to the case 
without time delay 0| (see Figs. [1] and [2]). In these 
figures, the value of A to reach the onset of synchroniza- 
tion varies periodically with r. For instance, in Fig. [H 
Ac(r = 1.0) < Ac(r = 0.3) < Ac(r = 1.5), while in Fig.l 
Ac(r = 1.0) < Ac(r = 0.3) < Ac(r = 2.0). 

Given this oscillatory dependence between O and r, we 
determine the value of r that makes the star networks 
more susceptible to the emergence of the onset of syn- 
chronization. More specifically, we found a time delay 
value that enhances the network synchronization level, 
while keeping the state stable in the thermodynamic li- 
mit K ^ oo. We also determined a value of time de- 
lay that yields the maximum value of critical coupling, 
i.e., a time delay that decreases the synchronization level. 
From equation Eq. [9l we obtained Ac = |cj — The- 
refore, while higher frequencies lead to a lower critical 
coupling enhancing the network synchronization, lower 
frequencies lead to higher ones, disturbing the network 
synchronization. In order to exam those extreme cases. 
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Let US now investigate more deeply the frequency de- 
pendency ^} = rt{r) supposing 



(11) 



where we have used the stationary solution with a phase 
difference Ohub,i = ^hub,it + (^hub^i- From Eq. [IH we 
can write Eq. [J] and [5] for the K -\- 1 nodes in the star 
network and thus leading to the following system of + 1 
equations 

K 

Q = K -\- X sin(— l^r + — an) 

ft = 1 — Xsm{^}r -\- ai — an) 
^} = 1 — X s'm{ftr + 0^2 — an) 



Figura 4: Coherence diagram r(A) for (a)K = 60, {h)K = 70, 
{c)K = 80, {d)K = 90. 




Figura 5: Comparison of the explosive synchronization obser- 
ved in star networks with different sizes K. The same delay 
shown in Figure [3] was adopted for the worst cases. The inset 
shows the same stars networks but for the case with r = 0. 



Figure 2] shows the coherence diagram for time delays 
chosen accordingly to Figure [3l For instance, for the 
time delay r = 2.4, the system presents the lowest possi- 
ble frequency, in such a way that the network undergoes a 
first-order transition to the synchronous state. Neverthe- 
less, for the time delay r = 0.2, no phase transition is 
observed and the order parameter r is increased until the 
synchronous state has been reached. 

Note that if the time delay r is set in order to establish 
a higher critical coupling Ac, a first-order phase transition 
is observed. Also, differently to the case in which r = 
12|, the dependence of the critical coupling A, 



with 



12| 





the number of nodes is lost. More specifically, in 
the authors show that the stability of the unlocked state 
(r ^ 0) increases as ^ oo, while here we observed that 
the transition to the synchronous state does not depend 
extensively on K (see Fig. [5]). 



ft = 1 — X sm{rtr + aK — (^h)- 

The average frequency Q of the star networks must satisfy 
all these equations. Summing all the equations in the 
system we have 

K 

{K^l)n = 2K^xY^ [s'm{-nr + a^) - sm{nr + a^)] . 

(12) 

Setting ttH = ai = a2 = • • • = On without loss of gene- 
rality we obtain 



2K 

-^[1-Asin(f]r)]. 



(13) 



According to [22|1 the system is linearly stable if and only 
if 

cos((7r) > (14) 

which means that the product must be in the range 

2n7r <Q.T < (2n + l)7r, n = 0, 1, 2, ... (15) 

Considering Q.T Qt — utt in Eq.[T3l in Figure[6]we show 
the solutions of Eq. [13] for different values of n. Observe 
that these solution correspond to the average frequency 
shown in Figure [3l Therefore, we verified a high agree- 
ment between theoretical and simulated results. Further- 
more, once we know exactly how the average frequency 
Q depends on r, is possible to determine the critical 
coupling for the star networks through the expression 
Ac = |cj — 1]| and with Eq. [T5]by rewriting Eq. [13] just in 
terms of Ac and r. Fig. [3 displays the results obtained for 
the critical couplings Ac for different values of r for the 
star network with K = SO, showing a perfect agreement 
for the theoretical values of Ac with the numerical data. 

In conclusion, the addition of time delay to the Ku- 
ramoto model allows reaching a higher synchronization 
level than those observed by Gardehes et al [l3|. Even 
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Figura 6: Theoretical and numerical curves for the frequency 
Q as function of r with A = 1.0, for the cases: (a)K = 60, 
(a)K = 70, {c)K = 80 and {d)K = 90. The discrete points 
correspond to the data shown in Figure [S] 



for a small time delay, the onset of synchronization is 
obtained for smaller coupling strengths than those for 
networks without time delay. The analysis performed 
here helps to understand real world communication sys- 
tems, in which the time delay is present due to the finite 
speed of the signal transmission across a communication 
medium. 
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